/**
 * @file   routinesO.cpp
 * @author Carlos Hernández García <carloshergar@usal.es>, Camilo Ruiz, Enrique Conejero
 * @date   Fri Apr  10 2009  
 *
 * Copyright (c) 2009 Carlos Hernández García
 *
 * @brief  See routinesO.h
 *
 * 
 * Carlos Hernández García                          \n
 * <carloshergar@usal.es>                           \n
 * Universidad de Salamanca                         \n
 * Plaza de la Merced s/n                           \n
 * E-37008 Salamanca, SPAIN                         \n
 */
/********************************************************************************/


#define c 137.037
#define pi 3.1416
#define nplim 2
#include <math.h>
#include <stdio.h>
#include <vector.h>
#include "../include/gridobject.h"
#include "../include/laserfield.h"
#include "../include/object.h"

/** 
 * This function calculates and distributes the charge of one particle into its nearest grid points
 * @version 1.0 (Enrique)
 * @param _part particle, whose charge is distributed
 * @param _cell pointer referencing the cell's grid
 * @param _qp array of charges of each particle
 * @param _deltag spatial grid increment
 * @param _na particle's size
 * @param _np number of particles
 * @param _ndg number of grid points
 */
void CalculaRo(particle *_part, cell1d *_cell,double _CalculaRofactor,double _deltag, int _na, int _np, int _ndg, FILE *fp0) //na is a parameter of the sheet
{
	double xa,xb,qa,qb,qi,qil;
	double xg;
	double qtotalfinal, qtotalinicial;
	
	/*qtotalinicial=0;
    for(int j=0;j<_xp.size();j++)
        qtotalinicial+=_qp[j];
    //printf("***Carga total antes de calculaRo= %.4f\n", qtotalinicial);*/
	
	
	for (int i=0;i<_ndg;i++)
    {
		_cell[i].ct.Rox=0.;
    }
	
	for(int j=0;j<_np;j++)
    {
		xg=_part[j].xp;
         
        int i=int(xg/_deltag);
		//fprintf(fp0, "%d, %e, %e, %e\n", i, xg, xg/_deltag, (xg/_deltag)+1);
		//printf("\nPosicion=%.3f\tPunto en la red (i)=%d\tCarga=%.4f\n", xg, i, _qp[j]);
        qa=_part[j].q*_CalculaRofactor*((i+1)*_deltag-xg)/_deltag;
		qb=_part[j].q*_CalculaRofactor*(xg-(i)*_deltag)/_deltag;
		//printf("qa=%.4f\tb=%.4f\n", qa, qb);
		
		xa=(xg+(i)*_deltag)/2.;
		xb=(xg+(i+1)*_deltag)/2.;
		qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
		
		for (int l=i-_na/2;l<=i+_na/2;l++)
		{
			if(l>=1 && (l<=_ndg))
			{	_cell[l].ct.Rox=_cell[l].ct.Rox+qi/_na/_deltag;
                //printf("qi:Asigno %.4f al punto %d\n", qi/_na/_deltag, l);
            }
		}
		qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
		//printf("qi=%.4f\tqi1=%.4f\n", qi, qil);
		for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
		{
			if( (l>=1) && (l<=_ndg))
			{	_cell[l].ct.Rox=_cell[l].ct.Rox+qil/_na/_deltag;
				//printf("qil:Asigno %.4f al punto %d\n", qil/_na/_deltag, l);
            }
		}
		
	}
	
}

/** 
 * This function calculates and distributes the current
 * @version 1.0 (Enrique)
 * @param _label flag indicating the axis direction (1:x, 2:y, 3:z)
 * @param _part particle, whose charge is distributed
 * @param _cell pointer referencing the cell's grid
 * @param _qp array of charges of each particle
 * @param _deltat temporal increment
 * @param _deltag spatial grid increment
 * @param _na particle's size
 * @param _np number of particles
 * @param _ndg number of grid points
 */
void CalculaJ(int _label, particle *_part, cell1d *_cell,double _CalculaRofactor,double _deltat,double _deltag, int _na, int _np, int _ndg, FILE *fp0, FILE *fp1)
{

	double xa,xb,qa,qb,qi,qil,q,mass;
	double xg;
	
	if(_label==1)
	{
		for(int i=0;i<_ndg;i++)
			_cell[i].vx.Jx=0.;
		
		for(int j=0;j<_np;j++)
		{
			xg=_part[j].xp;
			//int i=int(xg/_deltag)+1; MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double igamma=_part[j].igamma(); //sqrt(1.+(_part[j].vpx*_part[j].vpx+_part[j].vpy*_part[j].vpy+_part[j].vpz*_part[j].vpz)/lightC_au/lightC_au);
			
            //MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_part[j].q*_CalculaRofactor*((i+1)*_deltag-xg)/_deltag;
			qb=_part[j].q*_CalculaRofactor*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=1) && (l<=_ndg))
					_cell[l].vx.Jx=_cell[l].vx.Jx+qi/_na/_deltag*(_part[j].vpx)*igamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=1) && (l<=_ndg))
					_cell[l].vx.Jx=_cell[l].vx.Jx+qil/_na/_deltag*(_part[j].vpx)*igamma;
			}
			//fprintf(fp1, "%d\t%e\t%e\t%e\n",j, _vpx[j], _vpy[j], _vpz[j]);
			
		}//end j
		
    /*for (int i=0;i<_ndg;i++)
        fprintf(fp0, "[%d]: %e\n",i,_Jx[i]);*/
        
    //fprintf(fp0, "\n");		
       
		
	}//End lable1=Jx
	
	
	
	
	if(_label==2)
	{
		for(int i=0;i<_ndg;i++)
			_cell[i].ct.Jy=0.;
		
		for(int j=0;j<_np;j++)
		{
			xg=_part[j].xp;
			//int i=int(xg/_deltag)+1; //MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double igamma=_part[j].igamma();//double gamma=sqrt(1.+(_part[j].vpx*_part[j].vpx+_part[j].vpy*_part[j].vpy+_part[j].vpz*_part[j].vpz)/lightC_au/lightC_au);
			//MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_part[j].q*_CalculaRofactor*((i+1)*_deltag-xg)/_deltag;
			qb=_part[j].q*_CalculaRofactor*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=1) && (l<=_ndg))
					_cell[l].ct.Jy=_cell[l].ct.Jy+qi/_na/_deltag*_part[j].vpy*igamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=1) && (l<=_ndg))
					_cell[l].ct.Jy=_cell[l].ct.Jy+qil/_na/_deltag*(_part[j].vpy)*igamma;
			}
			
		}//end j
		
	}//end label2=Jy
	
	
	if(_label==3)
	{
		for(int i=0;i<_ndg;i++)
			_cell[i].ct.Jz=0.;
		
		for(int j=0;j<_np;j++)
		{
			xg=_part[j].xp;
			//int i=int(xg/_deltag)+1; //MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double igamma=_part[j].igamma();//double gamma=sqrt(1.+(_part[j].vpx*_part[j].vpx+_part[j].vpy*_part[j].vpy+_part[j].vpz*_part[j].vpz)/lightC_au/lightC_au);
			//MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_part[j].q*_CalculaRofactor*((i+1)*_deltag-xg)/_deltag;
			qb=_part[j].q*_CalculaRofactor*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=0) && (l<=_ndg))
					_cell[l].ct.Jz=_cell[l].ct.Jz+qi/_na/_deltag*(_part[j].vpz)*igamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=0) && (l<=_ndg))
					_cell[l].ct.Jz=_cell[l].ct.Jz+qil/_na/_deltag*(_part[j].vpz)*igamma;
			}
			
		}//end j
		
	}//end label==3
	
}//End CalculaJ

/** 
 * This function calculates tue incident electric field in the x direction
 * @version 1.0 (Enrique)
 * @param _xg position
 * @param _t time * @param _cell pointer referencing the cell's grid
 */
double E0x(double _xg, double _t)
{
	//E0x=0.01*sin(2.*w0*t)*sin(wo*t)*sin(w0*t);
	double ex=0.;
	return ex;
}

/** 
 * This function calculates tue incident gtric field in the y direction
 * @version 1.0 (Enrique)
 * @param _xg position
 * @param _t time * @param _cell pointer referencing the cell's grid
 * @param _xgleft initial position of the grid
 * @param _laser0 laser parameters
 */
double E0y(double _xg,double _t,double _xgleft, laser & _laser0, int _pulse_shape)
{
	double xx=_xgleft;
	double ey;
	
	if(_pulse_shape==0)
	{
		if((_xg-lightC_au*_t < xx) && (_laser0.kp*(_xg-xx)-_laser0.wp*_t >=-Pi) )
			ey=_laser0.E0*sin(_laser0.kp*(_xg-xx)-_laser0.wp*_t)*sin(_laser0.kp*(_xg-xx)-_laser0.wp*_t)*sin(_laser0.k0*(_xg-xx)-_laser0.w0*_t);
		else
			ey=0.;
	}
	return ey;
}

/** 
 * This function calculates tue incident electric field in the z direction
 * @version 1.0 (Enrique)
 * @param _xg position
 * @param _t time * @param _cell pointer referencing the cell's grid
 */
double E0z(double _xg, double _t)
{
	double ez=0.;
	return ez;
}

/** 
 * This function calculates tue incident magnetic field in the x direction
 * @version 1.0 (Enrique)
 * @param _xg position
 * @param _t time * @param _cell pointer referencing the cell's grid
 */
double B0x(double _xg,double _t)
{
	double bx=0.;
	return bx;
}

/** 
 * This function calculates tue incident magnetic field in the y direction
 * @version 1.0 (Enrique)
 * @param _xg position
 * @param _t time * @param _cell pointer referencing the cell's grid
 */
double B0y(double _xg,double _t)
{
	double by=0.;
	return by;
}

/** 
 * This function calculates tue incident magnetic field in the z direction
 * @version 1.0 (Enrique)
 * @param _xg position
 * @param _t time * @param _cell pointer referencing the cell's grid
 * @param _xgleft initial position of the grid
 * @param _laser0 laser parameters
 */
double B0z(double _xg,double _t,double _xgleft, laser _laser0, int _pulse_shape)
{
	double bz=E0y(_xg, _t, _xgleft, _laser0, _pulse_shape);
	return bz;
}


/** 
 * This function interpolates the electric and magnetic fields between the cells occuping the particle size (na)
 * @version 1.0 (Enrique)
 * @return electric and magnetic fields in the position of the particle
 * @param _x particle position
 * @param _cell pointer referencing the cell's grid
 * @param _Epx parameter that returns the electric field in the x-direction
 * @param _Epy parameter that returns the electric field in the y-direction
 * @param _Epz parameter that returns the electric field in the z-direction
 * @param _Bpx parameter that returns the magnetic field in the x-direction
 * @param _Bpy parameter that returns the magnetic field in the y-direction
 * @param _Bpz parameter that returns the magnetic field in the z-direction
 * @param _t time
 * @param _ndg number of grid points
 * @param _deltag spatial grid increment
 * @param _na particle's size
 */
void interpola(double _x,cell1d *_cell,double &_Epx,double &_Epy,double &_Epz,double &_Bpx,double &_Bpy,double &_Bpz, double _t,int _ndg,double _deltag, int _na, FILE *fp0)
{
	double xx=_x;
	//int i=int(xx/_deltag)+1;//MODIFICACION DE C
	int i=(int)round(xx/_deltag);
	_Epx=0.;
	_Epy=0.;
	_Epz=0.;
	_Bpx=0.;
	_Bpy=0.;
	_Bpz=0.;
	
	for(int l=i-_na/2;l<=i+_na/2;l++)
	{
		if((l>=0) && (l<_ndg) )
		{
			_Epx=_Epx+_cell[l].vx.Ex/_na;
			_Epy=_Epy+_cell[l].vx.Ey/_na;
			_Epz=_Epz+_cell[l].vx.Ez/_na;
			
			_Bpx=_Bpx+B0x(xx,_t)/_na;
			_Bpy=_Bpy+_cell[l].vx.By/_na;
			_Bpz=_Bpz+_cell[l].vx.Bz/_na;
		}
	}
	
}//end interpola



/** 
 * This function uses the Boris algorithm to move the particles, ie, since the fields are known,
 * calculates the velocity and the new position of the particles. (Based on the calculation of the relativistic Lorentz force)
 * @version 1.0 (Enrique)
 * @return particle velocity and position, inside the class particle called _part
 * @param _part particle
 * @param _Epx electric field in the x-direction acting in that particle
 * @param _Epy electric field in the y-direction acting in that particle
 * @param _Epz electric field in the z-direction acting in that particle
 * @param _Bpx magnetic field in the x-direction acting in that particle
 * @param _Bpy magnetic field in the y-direction acting in that particle
 * @param _Bpz magnetic field in the z-direction acting in that particle
 */
void ParticleMover(particle &_part, double _Epx,double _Epy,double _Epz,double _Bpx,double _Bpy,double _Bpz, double _deltat, double _lg, int &_coolparticles, int n, FILE *fp0)
{	
            double vvx,vvy,gvxs,gvys,vvz,gvzs,dcg,fb;
            double aexpt,aeypt,aezpt,abxpt,abypt,abzpt;
            double gamma;
			double ae=0.5*_part.qm*_deltat;
    
			//if(n>2800)fprintf(fp0,"\tMover entra: Epx=%f\tEpy=%f\tEpz=%f\t", _Epx, _Epy, _Epz);

			vvx=0.;
        	vvy=0.;
        	vvz=0.;
        	
        	gvxs=0.;
        	gvys=0.;
        	gvzs=0.;
        	
        	aexpt=0.;
        	aeypt=0.;
        	aezpt=0.;
			
			aexpt=ae*_Epx;
			aeypt=ae*_Epy;
			aezpt=ae*_Epz;
			
			abxpt=ae*_Bpx;
			abypt=ae*_Bpy;
			abzpt=ae*_Bpz;
			
			gvxs=_part.vpx+aexpt;
			gvys=_part.vpy+aeypt;
			gvzs=_part.vpz+aezpt;
			
			gamma=sqrt(1.+(gvxs*gvxs+gvys*gvys+gvzs*gvzs)/lightC_au/lightC_au);
			dcg=1./c/gamma;
			
			abxpt=dcg*abxpt;
			abypt=dcg*abypt;
			abzpt=dcg*abzpt;
			
			vvx=gvxs+gvys*abzpt-gvzs*abypt;
			vvy=gvys+gvzs*abxpt-gvxs*abzpt;
			vvz=gvzs+gvxs*abypt-gvys*abxpt;
			//if(n>2800)fprintf(fp0,"vvx=%f, vvy=%f, vvz=%f, gamma=%f\n", vvx, vvy, vvz, gamma);

			//Second half rotation and second half acceleration.
			fb=2./(1.+abxpt*abxpt+abypt*abypt+abzpt*abzpt);
			abxpt=abxpt*fb;
			abypt=abypt*fb;
			abzpt=abzpt*fb;
			
			gvxs=gvxs+vvy*abzpt-vvz*abypt+aexpt;
			gvys=gvys+vvz*abxpt-vvx*abzpt+aeypt;
			gvzs=gvzs+vvx*abypt-vvy*abxpt+aezpt;
			gamma=sqrt(1.+(gvxs*gvxs+gvys*gvys+gvzs*gvzs)/lightC_au/lightC_au);
			
			_part.vpx=gvxs;
			_part.vpy=gvys;
			_part.vpz=gvzs;
			
			_part.xpold=_part.xp;
			_part.xp=_part.xp+_deltat*_part.vpx/gamma;
			//Nos aseguramos de que nuestras particulas siguen dentro del rango en estudio (0,lg)
			if(_part.xp>=_lg)
			{	_part.xp=_part.xp-_lg;
				_coolparticles++;
				//printf("Particle out up!\n");
			}
			else if(_part.xp<0)
			{	_part.xp=_part.xp+_lg;
				_coolparticles++;
				//printf("Particle out down!\n");

			}	
			//if(n>2800)fprintf(fp0,"\tMover sale: xpold=%f\txp=%f\tvpx=%f\tvpy=%f\tvpz=%f\n", _part.xpold, _part.xp, _part.vpx, _part.vpy, _part.vpz);

}

/** 
 * This function calculates the internal fields due to the interaction between particles
 * @version 2.0 (Birdshall, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _deltat temporal increment
 * @param _xgleft init of the spatial grid
 * @param _laser0 laser parameters
 * @param _ndg number of grid points
 * @param _deltag spatial grid increment
 */
void CalculaCamposInternos(cell1d *_cell,double _deltat, double _xgleft, laser _laser0, int _ndg, double _deltag, FILE *fp0)
{
     
    //int ndg=_Ex.v.size();
    //double deltag=_Ex.dx1;
    double fp, fp_old, gm, gm_old, fm, fm_old, gp, gp_old;
    
    fp_old=_cell[0].vx.Fp;
    gm_old=_cell[0].vx.Gm;
    
    //Condiciones para la primera posicion del grid
    
    /***********************************************/
    //_Fp.v[0]= (E0y(0,_t,_xgleft, _laser0)+B0z(0,_t,_xgleft, _laser0))/2.;
    //_Gm.v[0]= (E0z(0,_t)+B0y(0,_t))/2.;
    
    _cell[0].vx.Fp= 0;
    _cell[0].vx.Gm= 0;    
    
    //cell->fp = laser_front.Qy * laser_front.field( time );
    //cell->gm = laser_front.Qz * laser_front.field( time + laser_front.shift );
    /***********************************************/                                   
    
    _cell[0].vx.Fm=_cell[1].vx.Fm- 2.*Pi*_deltat*_cell[1].ct.Jy;
    _cell[0].vx.Gp=_cell[1].vx.Gp- 2.*Pi*_deltat*_cell[1].ct.Jz;
    
	
    _cell[0].vx.Ey= _cell[0].vx.Fp+ _cell[0].vx.Fm;
	_cell[0].vx.Bz= _cell[0].vx.Fp- _cell[0].vx.Fm;
	_cell[0].vx.Ez= _cell[0].vx.Gp+ _cell[0].vx.Gm;
	_cell[0].vx.By= _cell[0].vx.Gp- _cell[0].vx.Gm;
    
	_cell[0].vx.Ex= _cell[0].vx.Ex-4.*Pi*_deltat* _cell[0].vx.Jx;
          
          
    for(int i=1; i<(_ndg-1); i++)
    {
          fp=_cell[i].vx.Fp;
          _cell[i].vx.Fp= fp_old- 2.*Pi*_deltat * _cell[i-1].ct.Jy;
          fp_old=fp;
          
          gm=_cell[i].vx.Gm;
          _cell[i].vx.Gm= gm_old- 2.*Pi*_deltat*_cell[i-1].ct.Jz;   
          gm_old=gm;                                       
          
          _cell[i].vx.Fm=_cell[i+1].vx.Fm- 2.*Pi*_deltat*_cell[i].ct.Jy;
          _cell[i].vx.Gp=_cell[i+1].vx.Gp- 2.*Pi*_deltat*_cell[i].ct.Jz;
    
          _cell[i].vx.Ey=_cell[i].vx.Fp+_cell[i].vx.Fm;
          _cell[i].vx.Bz=_cell[i].vx.Fp-_cell[i].vx.Fm;
          _cell[i].vx.Ez=_cell[i].vx.Gp+_cell[i].vx.Gm;
          _cell[i].vx.By=_cell[i].vx.Gp-_cell[i].vx.Gm;
          
          _cell[i].vx.Ex= _cell[i].vx.Ex-4.*Pi*_deltat*_cell[i].vx.Jx;
    }            
    
    
    //Condiciones para la ultima posicion del grid
    _cell[_ndg-1].vx.Fp= fp_old- 2.*Pi*_deltat*_cell[_ndg-2].ct.Jy;
    _cell[_ndg-1].vx.Gm= gm_old- 2.*Pi*_deltat*_cell[_ndg-2].ct.Jz;   
                                       
    /***********************************************/
    /*_Fm.v[ndg-1]= (E0y(ndg*_deltag,t,_xgleft, _laser0)-B0z(ndg*_deltag,t,_xgleft, _laser0))/2.;
    _Gp.v[ndg-1]= (E0z(ndg*_deltag,t)-B0y(ndg*_deltag,t))/2.;*/
    
    _cell[_ndg-1].vx.Fm= 0;
    _cell[_ndg-1].vx.Gp= 0;
    /***********************************************/
    
    _cell[_ndg-1].vx.Ey=_cell[_ndg-1].vx.Fp+_cell[_ndg-1].vx.Fm;
    _cell[_ndg-1].vx.Bz=_cell[_ndg-1].vx.Fp-_cell[_ndg-1].vx.Fm;
    _cell[_ndg-1].vx.Ez=_cell[_ndg-1].vx.Gp+_cell[_ndg-1].vx.Gm;
    _cell[_ndg-1].vx.By=_cell[_ndg-1].vx.Gp-_cell[_ndg-1].vx.Gm;
    
    _cell[_ndg-1].vx.Ex= _cell[_ndg-1].vx.Ex-4.*Pi*_deltat*_cell[_ndg-1].vx.Jx;
    
}


/** 
 * This function calculates the chrage deposition into the grid cells
 * @version 2.0 (Villasenor&Buneman, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _part particle whose charge we distribute into the cell's grid
 * @param _deltag spatial grid increment
 * @param _factor parameter to adecuate the lichters units to our atomic units, which is: nc*thickness/deltag/np
 */
void Charge_deposition (cell1d *_cell, particle &_part, double _deltag, double _factor)
{
	double here, next, prev;
	int i=int(_part.xp/_deltag);
	//double x0=(i+1)*_deltag;                   // contributions to charge of this cell, ...
	double x0=_cell[i].vx.x;
	double dist=(_part.xp - x0)/_deltag;
	
	//printf("part.xp=%f, i=%d, cell.vx.x=%f, dist=%f\n", _part.xp, i, _cell[i].vx.x, dist);
	//printf("x0(%d)=%f\tcell[%d].vx.x=%f\n", i, x0, i, _cell[i].vx.x);
	if (dist <= 0.5)
	{
		_cell[i].ct.Rox+= _part.zn*(0.5-dist)*_part.q*_factor;
		_cell[i-1].ct.Rox+= _part.zn*(0.5+dist)*_part.q*_factor;
	}	
	else 
	{
		_cell[i].ct.Rox+= _part.zn*(1.5-dist)*_part.q*_factor;
		_cell[i+1].ct.Rox+= _part.zn*(-0.5+dist)*_part.q*_factor;
	}
}



/** 
 * Function related to the calculation of the current deposition:
 * @version 2.0 (Villasenor&Buneman, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _part particle whose charge we distribute into the cell's grid
 * @param _deltag spatial grid increment
 * @param _factor parameter to adecuate the lichters units to our atomic units, which is: nc*thickness/deltag/np

 */
void left_one(cell1d *_cell,particle &_part, double _deltag,  double _factor,  int n, FILE *f)
{
	double xm=_part.xpold;               // before move
	double xp=_part.xp;                  // afterwards
	int i=int(xm/_deltag);
	double x0=_cell[i].vx.x;            // former left hand cell boundary

	//if(n>2800)fprintf(f, "\t\t\tleft_one, i=%d, entrada: Jx=%f\tJy=%f\tJz=%f\n", i, _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);
	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);

	
	double igamma=_part.igamma();
	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);

	//Cambios de unidades: multiplico por lightC_au y por "factor"
	double r_1  =0.5*_part.zn*(1.0 -(xp+xm-2.0*x0)/_deltag); 
	double r0   =0.5*_part.zn*(1.0 +(xp+xm-2.0*x0)/_deltag);
	
	//if(n>2800)fprintf(f, "\tigamma=%f\tr0=%f\t r_1=%f\t vpx=%f\tvpy=%f\tvpz=%f\n", igamma, r0, r_1, _part.vpx, _part.vpy, _part.vpz);
	

	
	_cell[i].vx.Jx+=(_part.zn*(xp-xm)/_deltag)*_part.q*_factor*lightC_au;
	_cell[i].ct.Jy+=r0*_part.vpy* igamma*_part.q*_factor;
	_cell[i].ct.Jz+=r0*_part.vpz* igamma*_part.q*_factor;

	if((i-1)>=0) 
	{
		_cell[i-1].ct.Jy+=r_1*_part.vpy* igamma*_part.q*_factor;
		_cell[i-1].ct.Jz+=r_1*_part.vpz* igamma*_part.q*_factor;
	}
	else
	{
		_cell[0].ct.Jy+=r_1*_part.vpy* igamma*_part.q*_factor;
		_cell[0].ct.Jz+=r_1*_part.vpz* igamma*_part.q*_factor;
	}
	//if(n>2800)fprintf(f, "salida: Jx=%f\tJy=%f\tJz=%f\n", _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);


}


/** 
 * Function related to the calculation of the current deposition:
 * @version 2.0 (Villasenor&Buneman, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _part particle whose charge we distribute into the cell's grid
 * @param _deltag spatial grid increment
 * @param _factor parameter to adecuate the lichters units to our atomic units, which is: nc*thickness/deltag/np
 */
void left_two_left(cell1d *_cell,particle &_part, double _deltag,  double _factor,  int n, FILE *f, int _ndg)
{
	double xm=_part.xpold;               // before move
	double xp=_part.xp;                  // afterwards
	int i=int(xm/_deltag);
	double x0=_cell[i].vx.x;                    // former left hand cell boundary
	
	//if(n>2800)fprintf(f, "\t\t\tleft_two_left, i=%d, entrada: Jx=%f\tJy=%f\tJz=%f\n", i, _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);
	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);
	double igamma=_part.igamma();

	//Cambios de unidades: multiplico por lightC_au y por "factor"
	double jx_1 = -_part.zn*(0.5 -(xp-x0+_deltag)/_deltag);
	double jx0  = -_part.zn*(0.5 +(xm-x0)/_deltag);

	double eps=(xm - (x0-0.5*_deltag) ) / ( xm - xp );
	double r_2=_part.zn* 0.5*(1.0-eps) * ( 0.5 - (xp-x0+_deltag)/_deltag );
	double r0=_part.zn* 0.5*eps * ( 0.5 + (xm-x0)/_deltag );
	double r_1=_part.zn - r0 - r_2;

	double jy_2 = r_2 * _part.vpy* igamma;
	double jy_1 = r_1 * _part.vpy* igamma;
	double jy0  = r0  * _part.vpy* igamma;
	double jz_2 = r_2 * _part.vpz* igamma;
	double jz_1 = r_1 * _part.vpz* igamma;
    double jz0  = r0  * _part.vpz* igamma;
 
	_cell[i].vx.Jx+= jx0*_part.q*_factor*lightC_au;
	_cell[i].ct.Jy+= jy0*_part.q*_factor;
	_cell[i].ct.Jz+= jz0*_part.q*_factor;
	if((i-1)>=0)
	{
		_cell[i-1].vx.Jx+= jx_1*_part.q*_factor*lightC_au;
		_cell[i-1].ct.Jy+= jy_1*_part.q*_factor;
		_cell[i-1].ct.Jz+= jz_1*_part.q*_factor;
	}
	else
	{
		_cell[_ndg-1].vx.Jx+= jx_1*_part.q*_factor*lightC_au;
		_cell[_ndg-1].ct.Jy+= jy_1*_part.q*_factor;
		_cell[_ndg-1].ct.Jz+= jz_1*_part.q*_factor;
	}
	if((i-2)>=0)
	{	_cell[i-2].ct.Jy+= jy_2*_part.q*_factor;
		_cell[i-2].ct.Jz+= jz_2*_part.q*_factor;
	}
	else
	{	_cell[_ndg-2].ct.Jy+= jy_2*_part.q*_factor;
		_cell[_ndg-2].ct.Jz+= jz_2*_part.q*_factor;
	}
	//if(n>2800)fprintf(f, "salida: Jx=%f\tJy=%f\tJz=%f\n", _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);

}

/** 
 * Function related to the calculation of the current deposition:
 * @version 2.0 (Villasenor&Buneman, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _part particle whose charge we distribute into the cell's grid
 * @param _deltag spatial grid increment
 * @param _factor parameter to adecuate the lichters units to our atomic units, which is: nc*thickness/deltag/np
 */
void left_two_right(cell1d *_cell,particle &_part, double _deltag,  double _factor,  int n, FILE *f, int _ndg)
{
	double xm=_part.xpold;               // before move
	double xp=_part.xp;                  // afterwards
	int i=int(xm/_deltag);
	double x0=_cell[i].vx.x;                    // former left hand cell boundary

	//if(n>2800)fprintf(f, "\t\t\tleft_two_right, i=%d, entrada: Jx=%f\tJy=%f\tJz=%f\n", i, _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);
	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);
	
	double igamma=_part.igamma();

    double jx0 = _part.zn * ( 0.5 - (xm-x0)/_deltag );
	double jx1 = _part.zn * ( 0.5 + (xp-x0-_deltag)/_deltag );

	double eps = ( x0 + 0.5*_deltag - xm ) / ( xp - xm );
	double r_1 = 0.5*eps * _part.zn * ( 0.5 - (xm-x0)/_deltag );
	double r1  = 0.5*(1.0-eps) * _part.zn * ( 0.5 + (xp-x0-_deltag)/_deltag );
	double r0  = _part.zn - r_1 - r1;

	double jy_1 = r_1 * _part.vpy* igamma;
	double jy0  = r0  * _part.vpy* igamma;
	double jy1  = r1  * _part.vpy* igamma;

	double jz_1 = r_1 * _part.vpz* igamma;
	double jz0  = r0  * _part.vpz* igamma;
	double jz1  = r1  * _part.vpz* igamma;

  
	if((i-1)>=0)
	{	_cell[i-1].ct.Jy+= jy_1*_part.q*_factor;
		_cell[i-1].ct.Jz+= jz_1*_part.q*_factor;
	}
	else
	{	_cell[_ndg-1].ct.Jy+= jy_1*_part.q*_factor;
		_cell[_ndg-1].ct.Jz+= jz_1*_part.q*_factor;
	}
	_cell[i].vx.Jx+= jx0*_part.q*_factor*lightC_au;
	_cell[i].ct.Jy+= jy0*_part.q*_factor;
	_cell[i].ct.Jz+= jz0*_part.q*_factor;
	if((i+1)<_ndg)
	{	_cell[i+1].vx.Jx+= jx1*_part.q*_factor*lightC_au;
		_cell[i+1].ct.Jy+= jy1*_part.q*_factor;
		_cell[i+1].ct.Jz+= jz1*_part.q*_factor;
	}
	else
	{	_cell[0].vx.Jx+= jx1*_part.q*_factor*lightC_au;
		_cell[0].ct.Jy+= jy1*_part.q*_factor;
		_cell[0].ct.Jz+= jz1*_part.q*_factor;
	}
	//if(n>2800)fprintf(f, "salida: Jx=%f\tJy=%f\tJz=%f\n", _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);

}

/** 
 * Function related to the calculation of the current deposition:
 * @version 2.0 (Villasenor&Buneman, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _part particle whose charge we distribute into the cell's grid
 * @param _deltag spatial grid increment
 * @param _factor parameter to adecuate the lichters units to our atomic units, which is: nc*thickness/deltag/np
 */
void right_one(cell1d *_cell,particle &_part, double _deltag,  double _factor,  int n, FILE *f, int _ndg)
{
  
	double xm=_part.xpold;               // before move
	double xp=_part.xp;                  // afterwards
	int i=int(xm/_deltag);
	double x0=_cell[i].vx.x;                    // former left hand cell boundary

	//if(n>2800)fprintf(f, "\t\t\tright_one, i=%d, entrada: Jx=%f\tJy=%f\tJz=%f\n", i, _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);
	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);

	double igamma=_part.igamma();

	double jx1 = _part.zn * (xp-xm)/_deltag;

	double r0  = 0.5 * _part.zn * ( 1.0 - (xp+xm-2.0*(x0+_deltag))/_deltag );
	double r1  = 0.5 * _part.zn * ( 1.0 + (xp+xm-2.0*(x0+_deltag))/_deltag );

    double jy0 = r0 * _part.vpy* igamma;
    double jy1 = r1 * _part.vpy* igamma;

    double jz0 = r0 * _part.vpz* igamma;
    double jz1 = r1 * _part.vpz* igamma;
    
	_cell[i].ct.Jy+= jy0*_part.q*_factor;
	_cell[i].ct.Jz+= jz0*_part.q*_factor;
	if((i+1)<_ndg)
	{	_cell[i+1].vx.Jx+= jx1*_part.q*_factor*lightC_au;
		_cell[i+1].ct.Jy+= jy1*_part.q*_factor;
		_cell[i+1].ct.Jz+= jz1*_part.q*_factor;
	}
	else
	{	_cell[0].vx.Jx+= jx1*_part.q*_factor*lightC_au;
		_cell[0].ct.Jy+= jy1*_part.q*_factor;
		_cell[0].ct.Jz+= jz1*_part.q*_factor;
	}
	//if(n>2800)fprintf(f, "salida: Jx=%f\tJy=%f\tJz=%f\n", _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);

}

/** 
 * Function related to the calculation of the current deposition:
 * @version 2.0 (Villasenor&Buneman, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _part particle whose charge we distribute into the cell's grid
 * @param _deltag spatial grid increment
 * @param _factor parameter to adecuate the lichters units to our atomic units, which is: nc*thickness/deltag/np
 */
void right_two_right(cell1d *_cell,particle &_part, double _deltag, double _factor,  int n, FILE *f, int _ndg)
{
	double xm=_part.xpold;               // before move
	double xp=_part.xp;                  // afterwards
	int i=int(xm/_deltag);
	double x0=_cell[i].vx.x;                    // former left hand cell boundary

	//if(n>2800)fprintf(f, "\t\t\tright_two_right, i=%d, entrada: Jx=%f\tJy=%f\tJz=%f\n", i, _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);
	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);

	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);
	
	double igamma=_part.igamma();

	double jx1 = _part.zn * ( 0.5 - (xm-x0-_deltag)/_deltag );
	double jx2 = _part.zn * ( 0.5 + (xp-x0-2.0*_deltag)/_deltag );

	double eps = (x0+1.5*_deltag - xm) / (xp - xm);
	double r0  = 0.5*eps * _part.zn  * ( 0.5 - (xm-x0-_deltag)/_deltag );
	double r2  = 0.5*(1.0-eps) * _part.zn  * ( 0.5 + (xp-x0-2.0*_deltag)/_deltag );
	double r1  = _part.zn - r0 - r2;
	
	double jy0 = r0 * _part.vpy* igamma;
	double jy1 = r1 * _part.vpy* igamma;
	double jy2 = r2 * _part.vpy* igamma;
	
	double jz0 = r0 * _part.vpz* igamma;
	double jz1 = r1 * _part.vpz* igamma;
	double jz2 = r2 * _part.vpz* igamma;
  

	_cell[i].ct.Jy+= jy0*_part.q*_factor;
	_cell[i].ct.Jy+= jz0*_part.q*_factor;
	if((i+1)<_ndg)
	{	_cell[i+1].vx.Jx+= jx1*_part.q*_factor*lightC_au;
		_cell[i+1].ct.Jy+= jy1*_part.q*_factor;
		_cell[i+1].ct.Jz+= jz1*_part.q*_factor;
	}
	else
	{
		_cell[0].vx.Jx+= jx1*_part.q*_factor*lightC_au;
		_cell[0].ct.Jy+= jy1*_part.q*_factor;
		_cell[0].ct.Jz+= jz1*_part.q*_factor;
	}
	if((i+2)<_ndg)
	{	_cell[i+2].vx.Jx+= jx2*_part.q*_factor*lightC_au;
		_cell[i+2].ct.Jy+= jy2*_part.q*_factor;
		_cell[i+2].ct.Jz+= jz2*_part.q*_factor;
	}
	else
	{
		_cell[1].vx.Jx+= jx2*_part.q*_factor*lightC_au;
		_cell[1].ct.Jy+= jy2*_part.q*_factor;
		_cell[1].ct.Jz+= jz2*_part.q*_factor;
	}
	//if(n>2800)fprintf(f, "salida: Jx=%f\tJy=%f\tJz=%f\n", _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);

}

/** 
 * Function related to the calculation of the current deposition:
 * @version 2.0 (Villasenor&Buneman, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _part particle whose charge we distribute into the cell's grid
 * @param _deltag spatial grid increment
 * @param _factor parameter to adecuate the lichters units to our atomic units, which is: nc*thickness/deltag/np
 */
void right_two_left(cell1d *_cell,particle &_part, double _deltag,  double _factor, int n, FILE *f, int _ndg)
{

	double xm=_part.xpold;               // before move
	double xp=_part.xp;                  // afterwards
	int i=int(xm/_deltag);
	double x0=_cell[i].vx.x;                    // former left hand cell boundary
  
	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);

	//if(n>2800)fprintf(f, "\t\t\tright_two_left, i=%d, entrada: Jx=%f\tJy=%f\tJz=%f\n", i, _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);
	//if(n>2800)fprintf(f, "\t\ti=%d xm=%f\n", i, xm);

	double igamma=_part.igamma();

	double jx0  = - _part.zn * ( 0.5 - (xp-x0)/_deltag );
	double jx1  = - _part.zn * ( 0.5 + (xm-x0-_deltag)/_deltag );

	double eps  = ( xm - (x0+0.5*_deltag) ) / ( xm - xp );
	double r_1  = 0.5*(1.0-eps) * _part.zn  * ( 0.5 - (xp-x0)/_deltag );
	double r1   = 0.5*eps * _part.zn  * ( 0.5 + (xm-x0-_deltag)/_deltag );
	double r0   = _part.zn  - r_1 - r1;

	double jy_1 = r_1 * _part.vpy* igamma;
	double jy0  = r0  * _part.vpy* igamma;
	double jy1  = r1  * _part.vpy* igamma;

	double jz_1 = r_1 * _part.vpz* igamma;
	double jz0  = r0  * _part.vpz* igamma;
	double jz1  = r1  * _part.vpz* igamma;
  

	if((i-1)>=0)
	{	_cell[i-1].ct.Jy+= jy_1*_part.q*_factor;
		_cell[i-1].ct.Jz+= jz_1*_part.q*_factor;
	}
	else
	{
		_cell[_ndg-1].ct.Jy+= jy_1*_part.q*_factor;
		_cell[_ndg-1].ct.Jz+= jz_1*_part.q*_factor;
	}
	_cell[i].vx.Jx+= jx0*_part.q*_factor*lightC_au;
	_cell[i].ct.Jy+= jy0*_part.q*_factor;
	_cell[i].ct.Jz+= jz0*_part.q*_factor;
	if((i+1)<_ndg)
	{	_cell[i+1].vx.Jx+= jx1*_part.q*_factor*lightC_au;
		_cell[i+1].ct.Jy+= jy1*_part.q*_factor;
		_cell[i+1].ct.Jz+= jz1*_part.q*_factor;
	}
	else
	{
		_cell[0].vx.Jx+= jx1*_part.q*_factor*lightC_au;
		_cell[0].ct.Jy+= jy1*_part.q*_factor;
		_cell[0].ct.Jz+= jz1*_part.q*_factor;
	}
	//if(n>2800)fprintf(f, "salida: Jx=%f\tJy=%f\tJz=%f\n", _cell[i].vx.Jx, _cell[i].ct.Jy, _cell[i].ct.Jz);


}

/** 
 * This function calculates the current deposition into the grid cells
 * We distinguish six cases:
 * First, distinguish former position in the first or second half of the cell
 * second, distinguish one one-boundary move and two two-boundary moves
 * #-boundary move means contributions to # boundary currents Jx.
 * The currents are calculated from the continuity equation,
 * assuming rectangular particle shape and area weighting
 * J.Villasenor and O.Buneman, Comp. Phys. Comm. 69 (1992) 306-316
 * @version 2.0 (Villasenor&Buneman, Lichters)
 * @param _cell pointer referencing the cell's grid
 * @param _part particle whose charge we distribute into the cell's grid
 * @param _deltag spatial grid increment
 * @param _factor parameter to adecuate the lichters units to our atomic units, which is: nc*thickness/deltag/np
 */
void Current_deposition(cell1d *_cell, particle &_part, double _deltag, double _factor, int n, FILE *f, int _ndg)
{

	double xm=_part.xpold;               // before move
	double xp=_part.xp;                  // afterwards
	int i=int(xm/_deltag);
	double x0=_cell[i].vx.x;                    // former left hand cell boundary
	
	double x0p05dx=x0 + 0.5*_deltag;
	double x0m05dx=x0 - 0.5*_deltag;
	double x0p15dx=x0 + 1.5*_deltag;
	
	//printf("curent_dep: x0=%f\txm=%f\tx0p05dx=%f\tx0m05dx=%f\tx0p15dx=%f\n", x0, xm, x0p05dx,x0m05dx,x0p15dx);
	
	if (xm<x0p05dx) 
	{                       // former position in first half of the cell
		if (xp<x0m05dx)    
		{	//if(n>21002)fprintf(f, "\t\t\t\tleft_two_left Jx=%f\tJy=%f\tJz=%f\n", i, _cell[(int)(_part.xpold/_deltag)].vx.Jx, _cell[(int)(_part.xpold/_deltag)].ct.Jy, _cell[(int)(_part.xpold/_deltag)].ct.Jz);
			left_two_left(_cell, _part, _deltag, _factor, n, f, _ndg);  // two boundary move to the left
		}
		else 
		{
			if (xp>=x0p05dx)
			{	//if(n>21002)fprintf(f, "\t\t\t\tleft_two_right Jx=%f\tJy=%f\tJz=%f\n", i, _cell[(int)(_part.xpold/_deltag)].vx.Jx, _cell[(int)(_part.xpold/_deltag)].ct.Jy, _cell[(int)(_part.xpold/_deltag)].ct.Jz);
				left_two_right(_cell, _part, _deltag, _factor, n, f, _ndg); // two-boundary move to the right
			}
			else                 
			{	//if(n>21002)fprintf(f, "\t\t\t\tleft_one Jx=%f\tJy=%f\tJz=%f\n", i, _cell[(int)(_part.xpold/_deltag)].vx.Jx, _cell[(int)(_part.xpold/_deltag)].ct.Jy, _cell[(int)(_part.xpold/_deltag)].ct.Jz);
				left_one (_cell, _part, _deltag, _factor, n, f);      // one boundary move
			}
		}
	}

	else 
	{                                 // former position in the second half of the cell
		if (xp>x0p15dx)    
		{	//if(n>21002)fprintf(f, "\t\t\t\tright_two_right Jx=%f\tJy=%f\tJz=%f\n", i, _cell[(int)(_part.xpold/_deltag)].vx.Jx, _cell[(int)(_part.xpold/_deltag)].ct.Jy, _cell[(int)(_part.xpold/_deltag)].ct.Jz);
			right_two_right(_cell, _part, _deltag, _factor, n, f, _ndg);// two boundary move to the right
		}
		else 
		{
			if (xp<=x0p05dx)
			{	//if(n>21002)fprintf(f, "\t\t\t\tright_two_left Jx=%f\tJy=%f\tJz=%f\n", i, _cell[(int)(_part.xpold/_deltag)].vx.Jx, _cell[(int)(_part.xpold/_deltag)].ct.Jy, _cell[(int)(_part.xpold/_deltag)].ct.Jz);
				right_two_left(_cell, _part, _deltag, _factor, n, f, _ndg); // two boundary move to the left
			}
			else                 
			{	//if(n>21002)fprintf(f, "\t\t\t\tright_one Jx=%f\tJy=%f\tJz=%f\n", i, _cell[(int)(_part.xpold/_deltag)].vx.Jx, _cell[(int)(_part.xpold/_deltag)].ct.Jy, _cell[(int)(_part.xpold/_deltag)].ct.Jz);
				right_one(_cell, _part, _deltag, _factor, n, f, _ndg);      // one boundary move
			}
		}
	}
}

double gauss_rand48(void)
{
  double r1, r2;

  r1 = drand48();
  r2 = drand48();

  return sqrt( -2.0 * log( 1.0 - r1 ) ) * sin( 2*Pi*r2 );
}


void InicioVelocidades(particle *_part, int VelocityInit, double Etherm, int _npion, double _vtherm, FILE *fp0)
{
    double pmax, dp;
    double pp, fv, df;
    double x[_npion];
    int iaux, j;
    double xaux;
    double vtherm;
    
    vtherm=(sqrt(2*Etherm))/lightC_au; //Etherm is given in u.a.
    printf("Etherm=%f\tvtherm=%f\n", Etherm, vtherm);
    
    if(VelocityInit==0) //Inicializamos todas las particulas en reposo
    {
        for (int j=0;j<_npion;j++)
        {
    		_part[j].vpx=0.;
    		_part[j].vpy=0.;
    		_part[j].vpz=0.;
	    }
	    return;
     }

     if(VelocityInit==1) //Inicializamos las  velocidades de las particulas segun BIRDSHALL    
     {              
        pmax=5.*sqrt(Etherm);
     	dp=2.0*pmax/_npion;
    	x[0]=0.0;
    	for (j=1;j<_npion;j++)
        {
    		pp=(j-0.5)*dp-pmax;
    		fv=exp(-pp*pp/(2.0*Etherm));
    		x[j]=x[j-1]+fmax(fv,0.0);
    	}
    	df=x[_npion-1]/_npion;
    	int i=0;
    	for (j=0;j<_npion;j++)
        {
    		fv=(+0.5+j)*df;
    		if(fv>=x[i+1]) 
            {
    			i++;
    		}
    		_part[j].vpx=dp*(i+(fv-x[i])/(x[i+1]-x[i]))-pmax;
    		_part[j].vpy=_part[j].vpx;
    		_part[j].vpz=0.;
    	}
    	for (j=0;j<_npion;j++)
        {
    		iaux=int(_npion*drand48());
    		xaux=x[j];
    		x[j]=x[iaux];
    		x[iaux]=xaux;
    		iaux=int(_npion*drand48());
    		xaux=_part[j].vpx;
    		_part[j].vpx=_part[iaux].vpx;
    		_part[iaux].vpx=xaux;
    		iaux=int(_npion*drand48());
    		xaux=_part[j].vpy;
    		_part[j].vpy=_part[iaux].vpy;
    		_part[iaux].vpy=xaux;
            //fprintf(fp0, "%d\t%e\t%e\n", j, _part[j].vpx, _part[j].vpy, _part[j].vpz);
    	}
    	return;
     }
     
     if(VelocityInit==2) //Inicializamos las  velocidades de las particulas segun LICHTERS    
     {
	   for (j=1;j<_npion;j++)
       {
           //do
		   //{
		      _part[j].vpx=_vtherm*gauss_rand48();//lightC_au;
		      _part[j].vpy=_vtherm*gauss_rand48();//lightC_au;
		      _part[j].vpz=_vtherm*gauss_rand48();//lightC_au;
		       //		  vx = exponential_rand( input.vtherm[part->species] ); vy = vz = 0.0;
	       //}
	       //while( _part[j].vpx*_part[j].vpx + _part[j].vpy*_part[j].vpy + _part[j].vpz*_part[j].vpz >= 1.0);         // make sure that |v| < 
           //fprintf(fp0, "%d\t%e\t%e\t%e\n", j, _part[j].vpx, _part[j].vpy, _part[j].vpz);
        }
        return;
     }  
}


/*
void Ionizacion(int _ndg, int _nplim,double _deltap0[],double _p0[],double _Ex[],double _Ey[],double _Ez[],double _wrate,double _t,double _deltat,int _ileft,int _iright,double _coef,double _exponente)//,int _ppp[][]
{
	double wrate[_ndg];
	for(int i=_ileft;i<=_iright;i++)
	{
		double aux=sqrt(_Ex[i]*_Ex[i]+_Ey[i]*_Ey[i]+_Ez[i]*_Ez[i]);
		if(aux>3.333e-3)
		{
			wrate[i]=_coef/aux*exp(-_exponente/aux);
		}
		else
			wrate[i]=0.;
		
		_deltap0[i]=_deltat*wrate[i]*_p0[i];
		if(_deltap0[i]>_p0[i] )
			_deltap0[i]=_p0[i];
		//_ppp[i][1]=0;
	}
}


void CalculaP0(double _xp[],double _vpx[],double _vpy[],double _vpz[],int _np,int _npion, int _ndg,double _qp[],int _nplim, double _ppp[][nplim+1],double _p0[],double _rop[],double _deltap0[], int n, int _iseed,double _deltag,int _ileft, int _iright,double _cargaminima,int _nppp,int _npminimo,int _npsinp)
{
	int jal=1;//+int(_npion*rand(_iseed));
	for(int j=jal;j<_npion;j++)
	{
		int i=int(_xp[j]/_deltag)+1;
		if( (i*_deltag-_xp[j])< (_xp[j]-(i-1)*_deltag )    )
			i=i+1;
		if(_ppp[i][1] < _nplim )
		{
			_ppp[i][1]=_ppp[i][1]+1;
			_ppp[i][1]=j;//+_ppp[i][1]]=j;  //Check this line
		}
	}
	for(int j=1;j<jal-1;j++)
	{
		int i=int(_xp[j]/_deltag)+1;
		if( (i*_deltag-_xp[j]) < (_xp[j]-(i-1)*_deltag)  )
		{
			i=i+1;
		}
		if(_ppp[i][1]<_nplim )
		{
			_ppp[i][1]=_ppp[i][1]+1;
			_ppp[i][1]=j;// +_ppp[i,1]]=j;  fix this !!!
		}
	}
	
	for(int i=_ileft;i<_iright;i++)
	{
		if(_deltap0[i]> _cargaminima)
			_nppp=_ppp[i][1];
		if(_nppp!=0 && _npion>_npminimo)
		{
			_p0[i]=_p0[i]-_deltap0[i];
			_rop[i]=_rop[i]+_deltap0[i];
			double aux=_deltap0[i]*_deltag/_nppp;
			for(int j=1;j<_nppp;j++)
			{
				int jb=_ppp[i][1+j];
				double aux2=_qp[jb]-aux;
				if(aux2!=0)
				{
					_vpx[jb]=_vpx[jb]*_qp[jb]/aux2;
					_vpy[jb]=_vpy[jb]*_qp[jb]/aux2;
					_vpz[jb]=_vpz[jb]*_qp[jb]/aux2;
					_qp[jb]=aux2;
				}
			}//end for
		}//end if
		else if(_npion > _np)
		{
			_npion=_npion+1;
			_qp[_npion]=-_deltap0[i]*_deltag;
			_vpx[_npion]=0.;
			_vpy[_npion]=0.;
			_vpx[_npion]=0.;
			_xp[_npion]=_deltag*(i-1.);
			_p0[i]=_p0[i]-_deltap0[i];
			_rop[i]=_rop[i]+_deltap0[i];
		}
		else
			_npsinp=_npsinp+1;
	}
	
}


double ran1(int _idum,int _iff)
{
	int m=714025;
	int ia=1366;
	int ic=150889;
	double rm=1./m;
	double temp;
	int ir[97];
	int iy;
	
	if(_idum > 0 || _iff==0)
	{
		_iff=1;
		//_idum=mod(ic-_idum,m);
		for(int j=1;j<97;j++)
		{
			//_idum=mod(ia*_idum+ic,m);
			ir[j]=_idum;
		}
		//_idum=imod[ia*_idum+ic,m];
		iy=_idum;
	}
	int j=1+(97*iy)/m;
	//if(j>97 || j<1)
	//	break;
	
	iy=ir[j];
	double tempran1=iy*rm;
	//temp=imod(ia*_idum+ic,m);
	ir[j]=_idum;
	return tempran1;
}
*/

/*
void realft(double data,int n,int isignn)
{
	double pi=4.*atan(1.);
	double theta=pi/n;
	double cl=0.5;
	
	if(isign==1)
	{
		c2=-0.5;
		four1(data,n,1);
	}
	else
	{
		c2=0.5;
		theta=-theta;
	}
	double wpr=-2.*sin(0.5*theta)*sin(0.5*theta);
	double wpi=sin(theta);
	double wr=1.+wpr;
	double wi=wpi;
	n2p3=2*n+3;
	for(int i=2;i<n/2;i++)
	{
		int il=2*i-1;
		int i2=il+1;
		int i3=n2p3-i2;
		int i4=i3+1;
		double wrs=wr;
		double wis=wi;
		
		double hlr=c1*(data[i1]+data[i3]);
		double hli=cl*(data[i2]-data[i4]);
		double h2r=-c2*(data[i2]+data[i4]);
		double h2i=c2*(data[i1]-data[i3]);
		
		data[i1]=h1r+wrs*h2r-wis*h2i;
		data[i2]=h1i+wrs*h2i+wis*h2r;
		data[i3]=h1r-wrs*h2r+wis*h2i;
		data[i4]=-h1i+wrs*h2i+wis*h2r;
		double wtemp=wr;
		
		double wr=wr*wpr-wi*wpi*wr;
		double wi=wi*wpr+wtemp*wpi+wi;
	}
	
	if(isign==1)
	{
		h1r=data[1];
		data[1]=h1r+data[2];
		data[2]=h1r-data[2];
	}
	else
	{
		h1r=data[1];
		data[1]=c1*(h1r+data[2]);
		data[2]=c1*(h1r-data[2]);
		four1(data,n,-1);
	}
}


void four1(double data,int nn, int isign)
{
	double pi=4.*atan(1.);
	int n=2*nn;
	int j=1;
	
	for(int i=0;i<n;i++)
	{
			if(j>i)
			{
				tempr=data[j];
				tempi=data[j+1];
				
				data[j]=data[i];
				data[j+1]=data[i+1];
				data[i]=tempr;
				data[i+1]=tempi;
			}
			m=n/2;
			if(m>=2 && j > m)
			{
				j=j-m;
				m=m/2;
			}
			j=j+m;
			mmax=2;
			if(n>mmax)
			{
				istep=2*mmax;
				theta=2.*pi/(isign*mmax);
				wpr=-2.*sin(0.5*theta)*sin(0.5*theta);
				wpi=sin(theta);
				wr=1.;
				wi=0.;
				for(int m=0;m<mmax;m++2)
				{
					for(int i=m;i<n;i++step)
					{
						j=i+mmax;
						tempr=wr*data[j]-wi*data[j+1];
						tempi=wr*data[j+1]-wi*data[j];
						
						data[j]=data[i]-tempr;
						data[j+1]=data[i+1]-tempi;
						data[i]=data[i]+tempr;
						data[i+1]=data[i+1]+tempi;
					}
					wtemp=wr;
					wr=wr*wprwi*wpi+wr;
					wi=wi*wpr+wtemp*wpi+wi;
				}
				mmax=istep;
			}
	}
}

*/
